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An efficient simulation method is presented for Brownian fiber suspensions, wliicli includes both 
uncrossability of the fibers and hydrodynamic interactions between the fibers mediated by a meso- 
scopic solvent. To conserve hydrodynamics, collisions between the fibers are treated such that 
momentum and energy are conserved locally. The choice of simulation parameters is rationalised 
on the basis of dimensionless numbers expressing the relative strength of different physical pro- 
cesses. The method is applied to suspensions of semifiexible fibers with a contour length equal to 
the persistence length, and a mesh size to contour length ratio ranging from 0.055 to 0.32. For such 
fibers the effects of hydrodynamic interactions are observable, but relatively small. The non-crossing 
constraint, on the other hand, is very important and leads to hindered displacements of the fibers, 
with an effective tube diameter in agreement with recent theoretical predictions. The simulation 
technique opens the way to study the effect of viscous effects and hydrodynamic interactions in mi- 
crorheology experiments where the response of an actively driven probe bead in a fiber suspension 
is measured. 

PACS numbers: 



I. INTRODUCTION 

The dynamics of rods and semifiexible fibers are 
strongly influenced by their mutual uncrossability. Ex- 
amples include carbon nanotubes [l| , fd- virus S Sjl, H, 
and biologically relevant polymers such as actin [aSBBl 
and tubulin To, . Already at surprisingly low concen- 
trations, uncrossability in such systems leads to a tem- 
porary and anisotropic "cage" or tube from which the 
rod or fiber can only escape through anisotropic motion 
(rcptation) 12] or through collective motion, as exem- 
plified by the collective reorientation observed in sheared 
concentrated rod suspensions [H, [H, ■ 

Besides the mutual uncrossability constraint, the dy- 
namics of rods and fibers are also infiuenced by Brownian 
forces (due to random collisions with solvent molecules) 
and hydrodynamic interactions (His) mediated by the 
solvent. The role of His in entangled suspensions of 
Brownian rigid rods and semifiexible fibers has remained, 
with a few exceptions, largely unexplored. This is caused 
by the difficulty of treating Brownian dynamics, hy- 
drodynamics, and entanglements within one theoreti- 
cal framework [l6| . His are dominant in the di- 
lute and onset of the semidilute regime. For exam- 
ple, the scaling of the relaxation times of the normal 
modes (Rouse modes) in an unentangled bead-spring 
chain changes from Tp oc (N/p)'^ for a chain without His 
to Tp PC (N/p)^^"^ for a chain with His (Zimm scaling) 
[T^ . Il7j . Here N is the number of beads and mode p mea- 
sures correlated motion on a length scale of {N/p) beads. 
Also the diffusion and segmental dynamics of dilute DNA 
molecules are controlled by hydrodynamic interactions 
[Tsl [Toj . On the other hand, it is believed that His are 
effectively screened in very concentrated suspensions and 
to a certain extent also in semidilute suspensions in equi- 
librium situations [1, [l^, . The onset of the semidi- 



lute regime already occurs at lower concentrations for 
rigid rods than for flexible chains of equal contour length 
jl2l |. This corresponds to a smaller dynamic correlation 
length in a semidilute suspension of rigid rods than in an 
equally concentrated suspension of flexible chains. In- 
deed, Pryamitsyn and Ganesan have shown that the ef- 
fects of His in semidilute and concentrated suspensions of 
completely rigid Brownian rods (with aspect ratio up to 
20) are secondary relative to the steric interactions [20j . 
A detailed analysis shows that His modify the diffusion 
parallel to the rod, in agreement with theories of hydro- 
dynamic screening [2ll . |22| . In all probability, the im- 
portance of His is decreasing with increasing chain stiff- 
ness and/or increasing concentration, but it is difficult 
to predict in general under which conditions His can be 
neglected. 

The need to consider His becomes particularly impor- 
tant when considering non-equilibrium situations. There 
are various applications where fibers are dragged along 
by flow or where the fibers generate flow because they 
are dragged by an external field. Examples include 
fiow through microchannels |23|. sedimentation or elec- 
trophoresis of fibers [23.[25l.l26||. and active microrheology 
27, 28, 29, 30]. In active microrheology a colloidal bead 
is embedded in a medium and driven by magnetic or 
optical forces. The force-displacement response is mea- 
sured with the goal to locally measure the rheology of 
the medium. In case of a medium consisting of a fiber 
network, it is important for the interpretation of these 
experiments to understand the hydrodynamic coupling 
between fluid flow generated by the probe bead on the 
one hand, and the fiber network on the other hand. The 
work presented here is part of a long-term effort to gen- 
erate this understanding. Coupling between fluid flow 
and fiber dynamics may be especially important when 
the probe bead is smaller than the mesh size of the net- 
work j^. Even for probe sizes in between the mesh size 
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and the fiber contour length interesting new mechanisms 
may be observed |3l| . 

Computer simulations in which His, entanglements, 
and Brownian motion are treated on an equal footing 
may help in gaining insight in the dynamics of Brow- 
nian fiber suspensions. First, let us focus on hydrody- 
namics. To rigorously include His in a simulation re- 
quires a decomposition of the mobility tensor, which is 
typically an 0{N^) operation [32], although with cer- 
tain approximations (expanding force distributions along 
rods in Legendre polynomials and retaining only lower or- 
der terms) this can be made more efficient [2J, [H, [sl] . 
Another approach is to explicitly include the solvent. 
The large gap in time- and lengthscales between the 
solvent molecules and colloidal sized particles has led 
to the development of mesoscopic simulation techniques 
which avoid the computationally costly explicit treat- 
ment of every solvent molecule. Important develop- 
ments in this area are Lattice Boltzmann (when ex- 
tensions to allow for thermal fluctuations are included) 
al, [si, HE HI, Dissipative Particle Dynamics (DPD) 
and Multi-Particle Collision Dynamics (MPCD) 
U4l, M |43|, [il |45l|46|, |47|, [isl, MM IsH Is^ . The 
latter, in its original implementation [40|, is also known 
as Stochastic Rotation Dynamics (SRD). All these meso- 
scopic simulation techniques account for correlated mo- 
tion of the solvent which leads to long-range hydrody- 
namic interactions. 

Second, let us focus on the entanglements. Most ex- 
isting methods implement non-crossing by resorting to 
explicit repulsive interactions. The dynamics of rela- 
tively short non-crossing rods may be modeled by means 
of forcefields with ellipsoidal or spherocylindrical geom- 
etry [sl, 113], whereas non-crossing rods or chains are 
often modeled by representing them as a string of rela- 
tively hard beads with bonds that are sufficiently strong 
to make the crossing of two such chains energetically un- 
favourable 0, m, . Although popular for its simple 
implementation, the latter approach has two disadvan- 
tages. Firstly, a large amount of beads is needed to rep- 
resent very long or very thin fibers or chains. Secondly, 
the use of hard excluded volume interaction potentials 
necessitates small time steps to accurately integrate the 
equations of motion. This makes the calculation of the 
dynamics of long thin rods and fibers computationally 
very costly. 

A few off-lattice methods exist that implement non- 
crossing chains without resorting to explicit repulsive in- 
teractions. Examples include a Brownian dynamics ac- 
ception/rejection scheme by Ramanathan and Morse [13] 
and the 'twentanglement' method of Padding and Briels 
m, [!§]. Both methods, however, are based on Brow- 
nian dynamics without His. This means that solvent- 
mediated interactions between the embedded chain seg- 
ments are ignored. Rather, the segments feel a certain 
friction with a fictitious static background fluid, as well 
as random forces. 

I will describe an efficient simulation algorithm for non- 



crossing fibers that includes hydrodynamic interactions. 
The method presented here relies on the SRD method to 
establish His between fiber or chain segments. In SRD 
a solvent is represented by Ng ideal particles of mass m. 
After propagating the particles for a time Stc, the sys- 
tem is partitioned into cubic cells of volume (with a 
random grid shift to conserve Galilean invariance [4l|). 
The velocities relative to the center of mass velocity of 
each separate cell are rotated over a fixed angle around 
a random axis. This procedure conserves mass, momen- 
tum, and energy and yields the correct hydrodynamic 
(Navier-Stokes) equations, including the effect of ther- 
mal noise [40|. The solvent particles only interact with 
each other through the rotation procedure, which can 
be viewed as a coarse graining of particle collisions over 
time and space. For this reason, the particles should 
not be interpreted as individual molecules but rather as 
a Navier-Stokes solver that naturally includes Brownian 
noise. The fiber or chain segments will be coupled to 
this hydrodynamic solvent by also taking part in the ro- 
tation procedure. With appropriately chosen simulation 
parameters (48j , such an approach leads to correct hy- 
drodynamic behaviour of polymeric chains, as shown re- 
cently by Winkler et al. [ITI]. From the point of view 
of the latter work, this paper is an extension of the hy- 
drodynamic method to also include uncrossability of the 
chains. 

This paper is organised as follows. A simple chain 
model is introduced in section [Til The non-crossing algo- 
rithm is described in detail in section [1111 The choice of 
simulation parameters is rationalised in section HVl and a 
validation and some results of the method are given in 
section [Vl Conclusions are given in section [VH 



II. CHAIN MODEL 

In this work a fiber or chain is represented by a string of 
vertices located at positions Ri {i — 1, . . . , Ny), with each 
vertex carrying a mass M. The non-crossing algorithm 
described in the next section is generally applicable to 
any model in which the interactions between connected 
vertices are described by potential energy terms. The 
model fiber or chain can achieve the right compressibility 
and bending stiffness by associating a bonding potential 
energy with each bond and an angular potential energy 
with each bend between two successive bonds. Specifi- 
cally, the potential energy of a bond + with length 
= JRi+i - Rj] is given by 
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(1) 



Here K is the elastic modulus of the fiber or chain and 
Iq is the equilibrium distance between each successive 
vertex. Two successive bonds (i — l,i) and + 1) 
with unit bondvectors Ui_i = (Rj — Ri_i)/i?j_i j and 
Ui = (Ri+i — Ri)/i?i,i+i make an angle 0i at vertex i. 
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with cos^i = Ui_i • Ui. The potential energy associated 
with this angle is given by 

/(ft,) = fcsr^(l-cos0O- (2) 

This particular form is chosen for the relative computa- 
tional ease of calculating cos 9i (rather than 9i). If ^ Ip 
the angles will typically be small and the angular poten- 
tial reduces to ^kBT{lp/lo)9f . Such a potential ensures 
that the persistence length of the fiber or chain is equal 
to Ip, as desired. Note that more realistic (non-linear) 
bond and angle potentials, as well as torsional stiffness 
effects, may be included but are ignored for simplicity. 

III. ALGORITHM 

In order to prevent chain crossing, a rule for the de- 
tection and treatment of bond-bond collisions needs to 
be devised. If hydrodynamic behaviour is to be retained, 
this collision rule must be compatible with the laws of 
conservation of momentum and energy. The most de- 
tailed approach would be to construct an event-driven 
algorithm in which a list of possible future collisions is 
generated and time progresses discretely from one colli- 
sion instant to the next. It is possible, but rather cum- 
bersome, to combine such a variable timestep algorithm 
with the fixed timestep SRD algorithm. However, re- 
solving the collisions to such detail is not in the same 
spirit as the SRD algorithm. In SRD one does not spec- 
ify the exact locations of the collisions between the sol- 
vent particles, but attains a rather more coarse grained 
view: collisions take place anywhere within the volume 
of a collision cell, anytime during the collision time in- 
terval. Technically, during the collision step the solvent 
particles are not actually displaced, they only exchange 
momentum and energy. This has proven to be sufficient 
for hydrodynamic behaviour of the solvent. 

In this work a similar fixed timestep idea is used for 
the collisions between chain bonds. It is unnecessary to 
specify the exact locations of the collisions. Rather, chain 
vertices are picked in random order and moved accord- 
ing to their velocities, except if this motion results in a 
collision with another chain [s^]. In the latter case mo- 
mentum and energy are exchanged between the vertices 
surrounding the colliding bonds. By moving the chain 
vertices one-by-one instead of all at once, the detection 
and treatment of the collisions are greatly simplified at 
the cost of accuracy in the collision location. It is neces- 
sary to use a random permutation for the order in which 
the vertices are picked, because otherwise bias may be 
introduced in successive collisions between the same pair 
of bonds. The solvent particles of mass m located at po- 
sitions Tj [j = 1, . . . , iVg) are treated as usual in SRD. 
Both solvent particles and chain vertices take part in the 
grid cell based collision step; this ensures that the chains 
are hydrodynamically coupled to the solvent. The algo- 
rithm may be summarised as follows: 



1. Read in coordinates (r,R) and velocities (v, V) of 
the solvent particles and chain vertices. 

2. Advance solvent positions over a timestep 5t 

Vj Yj + VjSt. (3) 

Apply periodic boundary or wall conditions to sol- 
vent coordinates. 

3. Create a randomly permuted list of all vertices. Try 
moving chain vertex from this list according to 

R, ^ Rt"al ^ + Y,6t- (4) 

Check for crossing of the bond {i—l,i) with another 
bond. Do the same for the bond {i,i + l). If a chain 
crossing occurs then reject this move, but exchange 
momentum and energy with the first collision part- 
ner. If no chain crossing occurs then accept this 
move. Apply periodic boundary or wall conditions 
to chain vertex coordinates. Details of crossing de- 
tection and momentum and energy exchange are 
given below. 

4. (May be performed less frequent:) The SRD col- 
lision step. Create a random-shifted grid and per- 
form random collisions of solvent and chain vertices 
within each grid cell according to 

*' 7 ' ^ * cm 

+ 7^(v, -V„„), (5) 

V, ^ Ycm+n{Y,-Vcm). (6) 

Here Vcm is the centre-of-mass velocity of all sol- 
vent and vertex particles in that particular cell and 
7?, is a rotation matrix which rotates velocities by 
a fixed angle a around a randomly oriented axis. 
Rescale velocities relative to centre-of-mass veloc- 
ity if thermostatting is required. 

5. Calculate vertex- vertex potential forces and possi- 
bly body forces for all particles: F.^ and ij. 

6. Advance velocities of solvent and vertices based on 
forces 

V, i-> V, + -^St (7) 
m 

V, ^ V, + ^St (8) 

7. If the number of required time steps has not yet 
been reached, go to step 2. 

8. Save coordinates and velocities of the solvent par- 
ticles and chain vertices. 

Most of the above algorithm is standard for SRD (note 
that in this version a leap-frog Verlet algorithm is used 
[iOl), except for step 3. If the update of the positions of 
the chain vertices would be treated similarly to step 2, 
then chains would be able to cross. More details on step 
3 are given in the next subsections. 
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FIG. 1: When vertex i is moved along its velocity vector 
Vi, a possible crossing of the connected bond + 1) with 
any of the neighbouring bonds is checked. The same applies 
to the connected bond (i — In this example a crossing 

between bonds + and + takes place. Momentum 
is exchanged along the direction n perpendicular to both these 
bonds at the time of impact. Note that in practice much 
smaller displacements of the vertices are used than shown 
here. This exaggerated view is only for reasons of clarity. 



A. Detecting bond crossings 

When performing a trial move of vertex i from R; to 
p^triai _ . _|_ Vi(5t, two bonds will move: {i — and 
+ 1) (see Fig. [1]). We assume that vertex i moves 
linearly in time, like 

R,(i) = R^ + V,^, te[o,6t]. (9) 

Focusing first on the bond + an intersection of this 
bond with another bond (j, j + 1) occurs at time tj if the 
vectors (Ri+i — Ri(t/)), (R^+i — R^) and (Rj+i — Rj) 
all lie within the same plane, i.e. if 

(Ri+i - R, - V,tj) ■ [(R,+i - Rj) X (R,+i - R,)] = 0. 

(10) 

This may be rewritten to 

^ (Rj+i - Ri) ■ [(Ri+i - Rj) X (Rj+i - Rj)] , , 
' V, • [(R,+i - Rj) X (R,+i - R,)] ■ ^ ' 

If tj e [0, 6t\ a collision may have occurred. Two further 
checks are needed to establish whether a real collision 
took place between the finite size bonds. If time is pro- 
gressed to the time of intersection tj, then points R(s) 
on bond {i, i + I) and R'(s') on bond + 1) are given 

by 

R(s) = R,(tj) + s(R,+i-R,(t,)), se[0,l](12) 
R'(s') = Rj + s' (Rj+i - R,) , s' e [0, 1] (13) 

The point of intersection, parametrised by the pair (s, s'), 
can be found by minimising the distance |R(s) — R'(s')| 



with respect to both parameters. The result is 



with 





s 


be — cd 
ac — b^^ 


(14) 
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ae — bd 
ac — b'^' 
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a — 


Rj+i - 




(16) 


b = 


(Rj+i ~ 


R-i(^/)) • ^ R-j) 7 


(17) 


c = 


R-j+i - 


Rjl^ , 


(18) 


d = 


(Ri+i - 


R,(t/))-(R,(t/)-Rj), 


(19) 


e = 




R,)-(R,(t/)-Rj). 


(20) 



Only if both s G [0, 1] and s' G [0, 1] a collision has oc- 
curred between the two finite bonds, and it occurred at 
time tj. 

A similar treatment is given to the bond pair (i — 1, i) 
and + 1). All neighbouring bonds + 1) which 
are not directly linked to the bonds {i — 1, i) or («, z -I- 1) 
must be checked in this way. The use of a Verlet linked 
list [13] greatly improves the efficiency of this procedure. 

If multiple collisions occur during the time interval 
[0, 6t\ due to the motion of a certain vertex i, the first 
collision is chosen for the exchange of momentum and 
energy, as discussed in the next subsection. The ratio of 
the number of executed collisions to the number of pos- 
sible collisions is monitored during the simulations. The 
integration time step should be so small that this ratio 
is close to one. 



B. Momentum and energy exchange 

Suppose that, as a consequence of the trial move of 
vertex i, a certain pair of bonds + 1) and + 1) 
have collided (the case of colliding bonds {i — and 
(ji j + 1) can be treated in a similar way). At the time of 
collision, tj, an amount AP of momentum is transferred 
from bond (i, z -|- 1) to bond (j, j -|- 1). This momentum 
transfer is directed along the normal to both bonds, i.e. 
AP = APn, with (see Fig. [T]) 

^ (R,+i-R,(tj))x (R,-+i-R,-) 

|(R,+i -R,(ij)) X (R,+i -R,)|- ^ ^ 

Note that in the simulation colliding bonds are not ac- 
tually moved (only non-colliding bonds are). The above 
calculation is needed to determine the direction in which 
momentum transfer is taking place. Because in this 
model the mass is concentrated in the vertices at the 
extremes of the bonds, the momentum transfer must be 
divided between the vertices following a lever rule. Us- 
ing the fact that all vertices have the same mass M, the 
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velocity change of the four vertices involved is given by: 

AV, = -(l-s)^n, (22) 
AP 

AV,+i 



AP 
— s n, 

M ' 

AP 

AV, = (l-,s')^n, 



M 



AV 



s n. 

M 



(23) 
(24) 
(25) 



Here s and s' are the fractional positions along the bonds 
where the collision has taken place. Note that this colli- 
sion automatically fulfills the law of conservation of mo- 
mentum. The amount of momentum transfer AP can 
be found from the law of conservation of energy. Before 
the momentum transfer the kinetic energy of the four 
involved vertices is given by 



(26) 



whereas after the collision it is given by 

2 



1 



aftc 



^ ' M 



V 



4 + 1 



AP, 

s n 

M 



AP 

V. + (l-.')^n 



Equating Xbcforo = ^aftcr we find 



(27) 



AP = 2M 



[jl - s)V, + sV,+i - (1 - s')Y, - s'Yj+,] ■ n 
(1 - s)2 + s-^ + {l- s')2 + s'2 

(28) 

to update the vertex 



This is used in Eqs 
velocities I6li. 



C. Extension to excluded volume fibers and chains 

The above method takes into account collisions be- 
tween infinitely thin fibers or chains. In some cases, for 
instance when the volume fraction is relatively large, it is 
desired to take into account the excluded volume of the 
fibers or chains. In this paper I will focus on the semidi- 
lute, low-volume fraction case where excluded volume is 
relatively unimportant (for example, the volume fraction 
will be such that no spontaneous nematic ordering will 
occur in the equivalent experimental system). However, 
for completeness, here follows an outline of the changes 
that need to be made to the algorithm; a detailed account 
will be given in a separate paper. 

When dealing with excluded volume it is envisaged 
that each bond + 1) represents the centre line of a 
tube of diameter D. The tube stretches from i to i -I- 1. 
Because the next bond {i + l,i + 2) is oriented differ- 
ently, one needs to be careful at the corners. This may 



be done by envisaging spheres of diameter D to be placed 
at the vertices. When moving vertex i, the detection of 
bond crossings is more complex than the case of thin 
lines because the time of collision cannot be determined 
independently from Eq. (jlOp anymore. Rather, a gener- 
alisation of Eq. is needed to indicate a point R(s; t) 
on the centre-line of bond («, i + I) at time t: 



R(s; t) = R, + Y,t + s(R,+i - R, - Y,t). 



(29) 



Eq. (fT5|) is still valid to indicate a point R'(s') on the 
centre-line of bond (j, j -I- 1) because this bond is not 
moved. Now multiple kinds of possible collisions need to 
be checked: between two bonds, between a bond and a 
vertex, and between vertices. The collision that has ac- 
tually taken place (if any within the indicated interval) is 
the one with the smallest associated collision time. These 
collision times are determined as follows. When checking 
bond + with (j, j + 1), the closest distance c?b6(t) is 
determined by functionally minimising |R(s;i) — R'(s')| 
with respect to the parameters s and s'. The time of 
impact then follows from d\jb{ti) — D. When checking 
bond + 1) with vertex j, the closest distance dtvit) 
is determined by functionally minimising |R(s;t) — Rj| 
with respect to the parameter s. The time of impact 
then follows from dbv{ti) = D. Finally, when checking 
vertex i with vertex j, the closest distance dy^{t) is given 
by dyy(t) = |Ri 4- Yit — Rj|. The time of impact then 
follows from dvv{ti) ~ D. Note that in all these cases a 
grazing collision could lead to two solutions of tj within 
the interval [0, <5<]. In that case the smallest of the two 
must be considered, as that will correspond to the incom- 
ing collision. 



IV. CHOICE OF PARAMETERS 

Before a system of semiflexible fibers or chains in a 
solvent can be simulated, a number of parameters need 
to be chosen. A summary of these parameters is given 
in Table HI In this paper lengths will be in units of cell 
size oq, energies in units of ksT , and masses in units of 
m (this corresponds to setting oq = 1, fc^T = 1, and 
TO = 1). Time, for example, is expressed in units of 

— oiQ \J m/kBT\ other units can be found in Table [H 
The exact values of the parameters will depend of course 
on the particular application in mind, but there are a few 
general rules which I will present here. 



Hydrodynamic coupling between the chains and 
the solvent 



The simulation method is supposed to capture the 
hydrodynamic interactions between different (parts of) 
chains. It is therefore important, first, to ensure that the 
solvent exhibits liquidlike momentum transfer, and sec- 
ond to ensure a sufficiently strong coupling between the 
chain vertices and the solvent. 
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TABLE I: Units and simulation parameters for semiflexible 
chains in an SRD fluid. The parameters listed in the table all 
need to be independently fixed to determine a simulation. 



Basic Units 
ao = length 
fcsT = energy 
m — mass 



Derived Units 



to — ao 



fcsT 



to 



Vo 



m \/mkBT 
taao al 



= time 

= diffusion constant 
viscosity 



Independent fluid simulation parameters 
7 = number of particles per cell 
5tc — SRD collision time step 
Q — SRD rotation angle 
Li, — box length 



Independent chain simulation parameters 



St = MD time step 

L = contour length 

lo = distance between successive vertices 

M = mass of a chain vertex 

K = elastic modulus 

Ip = persistence length 

Nc = number of chains 



Momentum transfer in the solvent is determined by 
the average number of fluid particles per cell 7, the time 
interval between collisions 6tc, and the rotation angle a. 
The simplicity of SRD collisions has facilitated the an- 
alytical calculation of many transport coefflcients of the 
solvent [13, m, HE 1131 ■ These analytical expressions are 
particularly useful because they enable an efficient tuning 
of the viscosity and other properties of the fluid, without 
the need for trial- and-error simulations. The viscosity 
has two contributions, kinetic and coUisional: 



Vcol 



jkBTStc 



57 



(7- H-e-T)(4- 2cosa- 2cos2a) 
m(l — cos a) 



ISaoStc 



-(7-1 + e-^) 



(30) 
(31) 



The kinetic viscosity must not be confused with the kine- 
matic viscosity i'. The latter, defined as v = rj/p = 
{rjkin + Vcol) / (m'y) , may be interpreted as the diffusion 
coefficient for momentum. In a liquid momentum diffu- 
sion is much faster than the self-diffusion of the solvent 
or solute molecules (the dimensionless Schmidt number is 
large [49l|). In SRD this may be ensured by choosing the 
collision time interval such that the mean free path be- 
tween collisions is at least one order of magnitude smaller 
than the collision cell size Oq, i.e. 6tc < O-Uq. For a de- 
tailed treatise the reader is referred to [49|. The tests 
described in the next section use 6tc = 0.02to- 

The vertices of the chain are coupled to the solvent 
by participating in the collision step. RipoU et al. (isj 
have shown that an optimal hydrodynamic coupling is 
achieved when the mass of the vertex is about equal to 
the total mass of the solvent particles in a cell and, as 
above, when the collision interval is chosen sufficiently 
small. The tests described in the next section use 7 = 5 
and M = 5m. Under these conditions the selfdiffusion 
Dm of the vertex is for a large part determined by hy- 
drodynamic correlations in the solvent. The effective 
hydrodynamic radius, deflned as ah = kBT/{6TTriDM), 
is approximately 0.3ao. The hydrodynamic interactions 
between different (segments of) fibers will then be cor- 
rectly reproduced if the equilibrium distance Iq between 
connected vertices is about twice the hydrodynamic ra- 
dius. Similar to the work of Winkler et al. we choose 
Iq = 0.5ao, which for flexible polymer chains was shown 
to yield the expected Zimm dynamics [17i |. 

The value of the rotation angle a also determines the 
amount of hydrodynamic coupling (isj . Obviously, the 
coupling will be less for smaller rotation angles; in the 
limit a = no momentum will be transfered between 
chain and solvent. Generally, in the range 7r/2 < a < tt 
the exact value of a is much less important for the cou- 
pling than the value of the collision interval (note that 
extremes near a = tt should be avoided). Since rota- 
tions around an angle of a — tt/2 can be implemented 
particularly efficiently, this value was chosen in all work 
described here. 



The SRD method has proven to be very robust when 
it comes to predicting hydrodynamic behaviour of em- 
bedded ob ject s, in both equilibrium and nonequilibrium 
situations S El, E El E 113, 111 lEl • The pre- 
cise speed of the dynamics depends on the choice of the 
above parameters, just as in a real experiment choosing 
glycerine instead of water will slow down the dynamics of 
embedded objects. Some choices will be computationally 
more efficient than others but as long as the appropri- 
ate limits mentioned above {Stc < O.lio, 7r/2 < a < tt, 
M w 7TO, and Iq « oq) are respected, the physical hydro- 
dynamic behaviour of the system will be correctly simu- 
lated. 



7 



TABLE II: Dimensionless numbers relevant to simulation of 
the statics and dynamics of a fiber suspension and the specific 
values used in the test simulations. When flow is applied (not 
in this work), one should be particularly mindful of the Mach 
and Reynolds numbers, which are reported here for the case 
of shear flow with shear rate 7. Recommended upper limits 
for Stokes flow are given. The Peclet number may be smaller 
or larger than 1, depending on experimental conditions. 



property 



definition 



value 



dimensionless persistence length Ip = Ip/L 

hydrodynamic aspect ratio p = L/{2ah) 

dimensionless mesh size ^* = L 

Compressibility effects Ma = "flp/cs 

Inertial vs. viscous forces Re = 7/^/1/ 

Convective vs. Brownian motion Pe = 7r,. 



B. Dimensionless numbers 



1 

64 

0.055 - 0.32 

< O.I 

< 0.1 



given by 



0.56(r 



(35) 



and a deflection length L^j (average distance between 
successive collisions of the fiber with its tube) Ld/L = 
0.64(r)^/^(^;)^/'^ + 0.39(r)*/^(^;)^/^- These expressions 
confirm the importance of the dimensionless numbers Z* 
and Note that Eq. (f35|) confirms the established seal- 



'1/5 



valid for long enough chains 



ing law L±^ 
65, 66, 67]. 

When flow is applied (this will be presented in a forth- 
coming article), a few more dimensionless numbers need 
to be taken into account to correctly characterise the 
relative strength of competing physical processes [11]. 
Firstly, the Mach number measures the ratio 



Ma: 



Vflow 



(36) 



Tuning of the model to experimental conditions is 
greatly facilitated by the use of dimensionless numbers. 
The dimensionless numbers which are relevant to a fiber 
suspension are summarised in Table [TTl The ratio of per- 
sistence length Ip to fiber contour length L, 



1* — In. 



(32) 



determines whether the fibers are fiexible (/* <C 1), semi- 
flexible {I* « 1) or stiff {Ip ^ 1). The aspect ratio of the 
flber 



L 

2ah ' 



(33) 



where ah is the (hydrodynamic) radius, is important for 
the hydrodynamic behaviour of the fiber. For example, 
the rotational and translational diffusion coefficient of a 
stiff rod depend strongly on p, even in dilute solutions [H, 
l62j . and the critical concentration for nematic ordering 
due to excluded volume depends on the ratio between 
persistence length and diameter, i.e. on pi* p^ . 

A network of fibers is further characterised by its mesh 
size: 



(34) 



Here c is the number density of fibers. The mesh size 
can be interpreted as an average distance between net- 
work segments, where the numerator 3 is a mere defini- 
tion. An important dimensionless number is the ratio of 
mesh size to contour length ^* = £,/L. Together with 
the dimensionless persistence length, it determines the 
amount of confinement that a fiber feels due to entan- 
glements with its neighbours [53, [H, HJ- For example, 
Hinsch et al. f64'| derive an effective tube diameter L 1 



the (relative) flow speed of the solvent. 

The 



between Vfi 

and c, — {b/'i)[kBT /m), the speed of sound 



Mach number measures compressibility effects since 
the sound speed is related to the compressibility of a 
liquid. It may sound obvious that Ma needs to remain 
small (<C 1) for physical fiber suspensions, but particle- 
based coarse-graining schemes drastically increase the 
Mach number. The fiuid particle mass m is typically 
much greater than the mass of a molecule of the un- 
derlying fluid, resulting in a lower speed of sound. In 
other words, particle based coarse-grained systems are 
typically much more compressible than the solvents they 
model. In practice, in order to avoid compressibility ef- 
fects in the dynamics of the system, the Mach number 
must remain lower than about 0.1 [49|. 

The Reynolds number is one of the most important di- 
mensionless numbers characterising hydrodynamic flows. 
Mathematically, it measures the relative importance of 
the non- linear terms in the Navier-Stokes equation [6^. 
Physically, it determines the relative importance of iner- 
tial over viscous forces and can be expressed as 



Re 



VflowR 



(37) 



where i? is a length scale relevant to the problem. For 
a flber suspension this could be the persistence length, 
i.e. R K, Ip. For micrometer sized objects, the Re is 
usually very small (Re <^ 1). The Reynolds number 
can be kept small by ensuring that the flow velocities 
do not exceed some maximum (this should be monitored 
during the simulation) and by choosing a relatively high 
kinematic viscosity. Again, the latter may be done by 
choosing a small collision interval Stc- 

The deflnitions of the Mach and Reynolds numbers 
above depend on the chosen relevant length scale as well 
as the characteristic flow velocity. In Table HIl we report 
Ma and Re for shear flow with shear rate 7, where the 
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relevant length scale of a semiflexible fiber is the persis- 
tence length, and the characteristic flow velocity is the 
maximal velocity difference over this length scale. 

Lastly, it is important that the relative importance of 
convective transport to diffusive transport is comparable 
between experiment and simulation. This is expressed 
by the Peclet number 



where D is the self-diffusion coefficient of the fiber. Alter- 
natively, under shear flow the Peclet number can be de- 
fined as the product of applied shear rate and the (Brow- 
nian) rotational relaxation time of the fiber or chain. 
In this respect it should be noted that the absence of ex- 
cluded volume interactions facilitates the simulation of 
very long and thin fibers, with very large characteristic 
times. For example, the characteristic times associated 
with rotational, perpendicular and parallel motion of a 
stiff rod scale like L^/[lnp + /(p)], with f{p) weak func- 
tions of the aspect ratio p [62]. In the large p limit, the 
friction perpendicular to the rod is twice that in the par- 
allel direction. This large p limit is reached (within a few 
percent accuracy) for p in the order of 30 [16[. Therefore, 
a connection between time in a simulation of rods with 
p = 30 and time in an experiment with much longer rods 
can be made by identifying the rotational relaxation time 
of the simulated rods with the rotational relaxation time 
in the experiment. 



C. Galilean invar iance 



the Mach number limit introduced above already lim- 
its the allowed flow velocities, and hence the magnitude 
of Vcm- Typically these limits imply |Vcm| < O.lao/^o 
and St < O.Olio, i.e. the error in the update of the 
centre-of-mass position of two colliding bonds is less than 
O.OOlflo. This is much smaller than any of the other typ- 
ical length scales (L, Ip, ^, ^q) of the problem. In the next 
section a test will show that the method is indeed effec- 
tively Galilean invariant for all tested flow velocities in 
the range < Vfiow < 0.54ao/to- 

D. Computational efficiency 

The ability to update the positions and check for col- 
lisions one vertex at a time makes the method efficient. 
Also, the use of MPCD to model the solvent makes the 
inclusion of hydrodynamic interactions relatively cheap. 
The precise speed of the simulation depends on the sys- 
tem size and chain density, where the computation rate 
scales approximately inversely linear with system volume 
and cL^. In its current implementation a system contain- 
ing about 1.6 X lO'^ solvent particles and 100 semiflexible 
fibers with an aspect ratio of p = 64 (i.e. represented by 
64 vertices each) at a density of cL^ — 100 is integrated 
at a rate of 30 time steps (St) per second on a modern 
single core processor. For this system one (dilute limit) 
rotational relaxation time r,. sa 8.3 x 10"* to is reached 
in 75 hours. Of course itself depends strongly on the 
length of the fiber. In the above example, when each 64 
vertex fiber is cut into two shorter fibers of 32 vertices, 
using the same mesh size, the time to reach decreases 
to 11 hours of computation. 



Although momentum is conserved locally in all sol- 
vent and bond collisions, the method presented here is 
not strictly Galilean invariant. Remember that a bond 
which exchanges momentum with its collision partner is 
not actually displaced. In this step reference is made to 
an absolute reference frame: the centre-of-mass of the 
collision partners should have been displaced over a dis- 
tance VcmSt, where Vcm is the centre-of-mass velocity 
of the collision partners. The influence of neglecting this 
centre-of-mass motion in one time step can be made ar- 
bitrarily small by choosing a sufficiently small molecular 
dynamics step St. When the number of time steps in 
which a particular bond collides is much smaller than 
the number of time steps in which the bond is moving 
according to its given velocity, correct dynamics is recov- 
ered. 

As it turns out, the above condition is not limiting 
the efficiency of the method, for three reasons. Firstly, 
the ratio of bond length to mesh size ^ is usually 
small, making collisions relatively rare for each partic- 
ular bond. Secondly, the molecular dynamics time step 
St already needs to be chosen relatively small to resolve 
the dynamics of the relative stiff springs needed to rep- 
resent real fibrillar materials, such as actin. Thirdly, 



V. VALIDATION AND RESULTS 
A. Dilute chains and fibers 

The dynamics of a flexible chain or semiflexible fiber 
in dilute solution is strongly affected by His. To test 
whether the SRD method indeed captures hydrodynamic 
interactions correctly, I will first focus on the qualita- 
tive and quantitative behaviour of the self-diffusion co- 
efficients of single chains or fibers [l3|- In all cases the 
solvent is represented by an average of 7 = 5 particles 
per cell. The collision interval is set to Stc = 0.02io and 
the collision angle to a — tt/2. 

Flexible chains are represented by = 5, 10, 20, 40, 80 
or 160 vertices of mass M = 5m at an equilibrium dis- 
tance Iq = and a bond strength k — ZksT j af^, cor- 
responding to an entropic spring with root-mean-square 
bond length / = lao |12l | (all angular interactions have 
been disabled). The size of the cubic periodic simulation 
box is varied linearly with the root-mean-square end-to- 
end distance Re — 1\/N to avoid artifacts due to finite 
system sizes. Explicitly, Lb = 25ao is chosen for N = 10. 
For flexible chains, hydrodynamic Zimm theory predicts 
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FIG. 2; Simulation results for the center of mass diffusion 
coefficient (in units ao/fo) for dilute flexible chains (black cir- 
cles) and rigid-like rods (Lp = 2L; red squares) for various 
contour lengths. The solid and dashed lines correspond to 
theoretical predictions Eqs. I|39|) and (|40p . The dot-dashed 
line indicates the scaling expected for chains or rods without 
hydrodynamic interactions. 



a self-diffusion coefficient given by [1 

.ksT 



D = 0.196 



(flexible chain) 



(39) 



Fig. [2] presents the diffusion coefficients of the centres of 
mass of flexible chains of various length (black circles). 
Qualitatively, a scaling D oc (x 7V~^/^ can be ob- 
served. Quantitatively, using the analytically known vis- 
cosity from Eqs. ((30)) and ([3T|) . good agreement is found 
if the prefactor 0.196 in Eq. ([55]) is replaced by 0.17 (solid 
line). A slightly lower self-diffusion is in agreement with 
the fact that periodic images of the chain interact with 
each other via the periodic boundaries [49| . 

For rigid rods, the self-diffusion coefficient is given by 

M 



D = 







Kf) 


+ 0.312 







(rigid rod) (40) 



(up to order b/L), where b is the hydrodynamic diam- 
eter of the rod. To verify this relation, single rod-like 
fibers are represented by = 10, 15, 30, 45 or 60 vertices 
of mass M = 5m at an equilibrium distance Iq = 0.5ao 
and bond springs with strength k = K/Iq ~ lOOfesT/oQ. 
In order to minimise effects of fiexibility, the persis- 
tence length is chosen equal to twice the contour length, 
Ip — 2L. The relatively stiff bonds and angles require a 
molecular dynamics integration step of 6t — O.Ol^o- To 
avoid artifacts in the determination of the fiber length 
dependence due to finite system size effects, the size of 
the cubic periodic simulation box is increased linearly 



with the length of the fiber, where Lb = ISoq is chosen 
for N = 10. Fig. [5] presents the diffusion coefficients of 
the centres of mass of rod-like fibers of various length 
(red squares), together with the theoretical curve Eq. HOI 
(dashed line). Similar to the work described in p/7| the di- 
ameter b and the prefactor are obtained by a least squares 
fit, yielding b = O.Gao and a prefactor 0.094. The diame- 
ter is in good agreement with the effective hydrodynamic 
radius estimated for our vertices. The prefactor is slightly 
smaller than the theoretical prediction l/(37r) = 0.106, 
which can again be attributed to the slowing effect of 
periodic images. 

Note that in the absence of hydrodynamic interac- 
tions each vertex would act as an independent source of 
friction, leading to a centre-of-mass diffusion coefficient 
which scales like (dot-dashed line) for both flexi- 

ble chains and rigid rods. From these tests it may be 
concluded that the SRD method correctly captures the 
hydrodynamic interactions for flexible chains and semi- 
flexible flbers. 



B. Semidilute fibers 

In the following tests I will focus on the dynamics of 
suspensions of many semiflexible fibers, each similar to 
the rod-like fibers studied above, but now represented by 
64 vertices and a persistence length equal to the contour 
length, i.e. Ip = L = 32ao. All simulations were per- 
formed in a periodic cubic box with sides Lf, = 32aQ. 
The number density c of fibers was varied between the 
values cL^ — 30, 100, 300 and 1000, corresponding to 
mesh sizes — 10.2ao, 5.44ao, 3.20ao and 1.76ao, respec- 
tively. Higher values of the fiber density are not relevant 
because excluded volume effects can then no longer be 
neglected [l^ . 



1. Validation of Galilean invariance 

To test the effective Galilean invariance of the non- 
crossing constraint explicitly, a periodic system of semi- 
flexible fibers at the highest density of cL^ = 1000 
was subjected to a homogeneous flow in the x-direction 
with velocities ranging from zero to a relatively high 
Vfiow = 0.54ao/io- During a run of 10^ integration steps, 
several properties were monitored and compared to a sys- 
tem at rest (vfiow = 0). 

The energy and the centre-of-mass velocity of the sys- 
tem was observed to remain exactly constant. This con- 
flrms that energy and momentum are conserved during 
the fiber collisions also in the presence of background 
flow. 

The vertex mean square displacement 



g(t) = ^(R,(t)-R,(0))^) (41) 
averaged over all vertices i, as well as the mean square 
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FIG. 3: Mean square displacement of the vertices g{t) and 
the fiber centres-of-mass Qcmit), normalised by the square of 
the contour length L of the fiber, measured relative to a back- 
ground fiowing with velocity Vfiow (in units ao/to) in the x- 
direction, in a system with Ip — L and cL^ — 1000. Time 
is normalised by the rotation time of a fiber in the dilute 
limit. The results are indistinguishable for all flow velocities 
up to 0.54ao/io, signifying effective Galilean invariance of the 
method. 



FIG. 4: Mean square displacement (ci'^(t)) of the minimum 
distance between the middle vertex at time t + r and the flber 
at time t, without the non-crossing constraint (unentangled, 
dashed lines) and with the non-crossing constraint (entangled, 
solid lines), and with hydrodynamic interactions (HI, black 
colored lines) and without hydrodynamic interactions (no HI, 
grey lines). The dash-dotted line shows the tube diameter 
predicted by Eq. (I35|) . For all data shown Ip = L and ci^ — 
1000. 



displacement 

(t) = ((Rc™(t)-Rc™(0))') (42) 

of the fiber centres-of-mass (both relative to the back- 
ground flow) were determined and observed to be nearly 
indistinguishable, as shown in Fig. [3l This conclusively 
shows that, for the chosen parameters, the method is ef- 
fectively Galilean invariant for all relevant flow velocities. 

2. Influence of hydrodynamic interactions and 
uncross ability of fibers 

Hydrodynamic interactions may easily be turned off by 
selecting random pairs of fluid particles after the collision 
step and exchanging their velocities. In this manner en- 
ergy and momentum are still conserved globally, but no 
longer locally. The non-crossing constraint can be turned 
off by simply skipping the bond collision check. 

Figure m shows the effect of hydrodynamic interactions 
and uncrossability of fibers on a quantity {(P{t)^, where 
d{t) is defined as 

rf(i) =min|R,„(t + T)-Rj(T)|, (43) 

3 

where m = is the middle vertex of a fiber and j 

runs over all vertices 1,...,A^„ of that fiber. In other 
words, d is the closest distance between the position R,„ 



of the middle vertex at time t-\-T and any of the vertices 
of the fiber at an earlier time t. In a tightly entangled 
solution, the magnitude of the plateau in this quantity 
is a measure of the width of the tube to which the fiber 
is confined [s^, [1^ . The time axis is normalised by the 
rotation time of a fiber in the dilute limit, measured 
from the end-to-end vector decorrelation of a single fiber 
in a box of the same dimensions and with or without His, 
respectively. Two observations can be made. 

First, the results without His (grey lines) are system- 
atically below the results with His (black lines). The rel- 
ative difference is larger at shorter correlation times than 
at longer correlation times, leading to small differences in 
scaling of Id^lt)) with time t. It may be concluded that, 
apart from such small differences, the overall behaviour 
with or without His is quite similar for semiflexible fibers 
of length L = Ip. This result is in agreement with find- 
ings for completely rigid rods (Ip S> L) where it was found 
that the effects of His are secondary relative to the steric 
interactions [20| . 

Second, the results using the non-crossing constraint 
(solid lines) are equal to the results without this con- 
straint (dashed lines) at short times, whereas they devi- 
ate significantly at larger times. The transition between 
these two regimes may be interpreted in the tube model 
(T^ as the moment when the fibers start to collide with 
their effective tube walls. Fig.|3]shows the prediction 2L\ 
of Eq. ([55)1 (horizontal dash-dotted line labeled 'tube'), 
where the factor of 2 arises because in the theory of Ref. 
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FIG. 5: Mean square displacement of the minimum 

distance between the middle vertex at time t + r and the fiber 
at time r with the non-crossing constraint and hydrodynamic 
interactions for concentrations cJJ" = 30, 100, 300 and 1000. 
The dashed line indicates the expected t^''* behaviour for a 
single semiflexible fiber. For all data shown — L. 

[S^ l L± is defined as the mean square transverse dis- 
placement of one Cartesian component only. The agree- 
ment between the observed and predicted tube diameter 
is good. 

Focusing now on the most realistic case, with His and 
non-crossing fibers, the influence of network density is 
shown in Fig. [Sj For the lowest density shown, cL^ = 30, 
the fibers behave almost as dilute single fibers. In this 
limit, the growth of with time t at early times is 

limited by the finite transversal fluctuations of a wormlike 
chain, leading to an expected i^/"* scaling [t^. This is 
indeed observed in the simulations as well (dashed line). 
With increasing network density (and hence decreasing 
mesh size) the displacement of the fibers become hindered 
by the presence of other fibers at smaller and smaller 
length scales. A more in-depth analysis will be presented 
in a forthcoming paper. 



VI. CONCLUSIONS 

I have introduced a method to simulate the dynam- 
ics of Brownian fiber suspensions, where hydrodynamic 
interactions are mediated by a mesoscopic solvent and 
collisions between fibers are treated such that momen- 
tum and energy are conserved locally. The method is 
made efficient by moving one fiber segment at a time in- 
stead of all segments at once. A similar idea was used 
in the work of Ramanathan and Morse 5f\ in the con- 
text of non-hydrodynamic Brownian dynamics, whereas 
in this work hydrodynamics are conserved. The effective 
Galilean invariance of the current method was explicitly 
checked. 

It was found that for semidilute semiflexible fibers with 
L — Ip the effects of hydrodynamic interactions are small 
compared to the effects of uncrossability of the fibers. 
Because a similar observation has already been made for 
completely rigid rods [l^l, it may be concluded that His 
are relatively unimportant for all semidilute suspensions 
of fibers for which L < Ip. This is also the reason why 
the observed displacements of fibers in a hydrodynamic 
solvent are globally similar to those obtained in non- 
hydrodynamic simulations [13, , although differences 
are observed upon closer inspection. At constant chain 
concentration, these differences will become increasingly 
more important for longer chains {L > Ip) (iTj . or in sit- 
uations where fibers are subjected to fiow. The purpose 
of this paper was to introduce and validate the method; 
in a forthcoming paper I will focus on non-equilibrium 
situations. For example, the effect of viscous drag and 
hydrodynamic interactions will be studied in microrhcol- 
ogy experiments where the response of an actively driven 
probe bead in a fiber suspension is measured. 
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